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We employ molecular dynamic simulations to study the reduction process of graphene-oxide (GO) 
in a chemically active environment enriched with hydrogen. We examine the concentration and 
pressure of hydrogen gas as a function of temperature in which abstraction of oxygen is possible 
with minimum damage to C-sp^ bonds hence preserving the integrity of the graphene sheet. Through 
these studies we find chemical pathways that demonstrate beneficiary mechanisms for the quality 
of graphene including formation of water as well as suppression of carbonyl pair holes in favor of 
hydroxyl and epoxy formation facilitated by hydrogen gas in the environment. 



I. INTRODUCTION 

Following the progress in the fabrication and me- 
chanical exfoliation of the graphene^ ^, extraordinary 
low-dimensional electrical and mechanical properties of 
graphene have been revealed^. Currently, mass-scale 
fabrication of graphene-based devices free from defects 
and imperfections requires chemical processing with high 
quality product and is one of the challenging technologi- 
cal problems for graphene device technology. For exam- 
ple, the micro- mechanical cleavage of graphite^, epitaxial 
growth on silicon carbide^, chemical vapor deposition of 
hydrocarbons on transition metal surfaces^ have shown 
difficulties in obtaining processable graphene sheets in 
large quantities thus impeding full exploitation of its ex- 
citing properties. Chemical oxidation of graphite and 
subsequent exfoliation in solution allow a large scale pro- 
duction of isolated graphene oxide (GO). However, the 
reduction of GO is known to create many structural de- 
fects degrading the material properties of reduced GO. 
Hence identifying reliable methods that allow removal of 
oxygen from GO with minimum damage to C-sp^ bonds 
in graphene-structures is the focus of current research 
activities. 

Recently studies have shown that the chemical reduc- 
tion of graphite oxide can be a low-cost and scalable 
methocP^^. It has been found that layered materials 
constituting graphene layers functionalized with epoxy 
and hydroxyl groups are easily exfoliated in water. The 
resulting graphene oxide (GO) monolayers can be de- 
posited in controllable density onto a large variety of sub- 
strates, thus enabling the preparation of thin conductive 
films on solid and flexible substrateJ^^H^. 

It is accepted that GO can be described as a ran- 
dom distribution of oxidized areas with the oxygenated 
functional groups, combined with nonoxidized regions 
wherein most of the carbon atoms preserve sp^ hybridiza- 
tioiP^. GO is electrically insulating material, however, 
its conductivity can improve up to 4 orders of magni- 
tude by chemical reductiorP^"^^, nevertheless the typical 
conductivities of reduced GO (RGO) is stil l lower than 
pristine graphene by a factor of 10-lOC^^' due to pres- 
ence of residual functional groups remaining after reduc- 



tion. In contrast to mechanically exfoliated graphene, the 
chemically derived graphene (i.e., RGO) is found to con- 
tain a considerable amount of topological defects. Sim- 
ilar to other low-dimensional carbon nanostructures like 
carbon nanotubes^^, graphene ribbonJ^^^^, fullereneP^ 
and more recently graphene quantum dots^^, topological 
defects and etch-holes (that are unavoidable products of 
the reduction process) are expected to strongly affect the 
device electronic and mechanical performance, and thus 
to account for the differences between RGO and pristine 
graphene. 



Recent advances in computational modeling of mate- 
rials such as molecular dynamic (MD) simulations and 
ab-initio calculation have enabled researchers to provide 
detailed studies in the microscopic level and allows a care- 
ful analysi s of th e chemical reactions and diffusion mech- 
anisms^^ ^-*^ ! ^^ ! ^^ !. Such studies are helpful in identifying 
the experimental conditions that allows the minimization 
of mechanical and chemical damages to graphene and the 
optimization of their performance^^ 



In this work we employ MD sim ulations based on ab- 
initio CPMEPi and ReaxFpE^ which are the mathe- 
matical formulation that governs the appropriate dynam- 
ics of the molecular system to analyze the carbon-oxygen 
chemical reaction and reduction process. We show that 
the GO structural relaxation assisted with the oxygen 
surface diffusion that allows the migration of epoxy func- 
tional groups and protects the C-sp^ bonds, preserve the 
integrity of the graphene lattice structure and governs 
the dynamical pathway in low temperatures. In con- 
trast, in high temperatures, consistent with recent find- 
ing of Ref.^^ we show that the formation of carbonyl 
is dynamically favorable and can be accounted for domi- 
nant mechanism in breaking C-sp^ bonds. We then con- 
sider immersed systems of GO in dilute gas of hydro- 
gen using two large-scale structural models for GO and 
demonstrate the effect of hydrogen in removal of oxygen, 
blocking the carbonyl pair formation and finally seek for 
conditions in which structural damage to C-sp^ bonds is 
minimized. 
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II. MODEL AND COMPUTATIONAL 
METHODS 

The atomic model structure of large scale GO is con- 
structed by including randomly distributed hydroxyl and 
epoxy functional groups to a single-graphene sheet. The 
dynamical stability of the entire structure in particular 
the functional groups attached to the graphene will be in- 
vestigated by performing MD in which the time evolution 
of the system and its environment is calculated. The sys- 
tem consists of GO in vacuum or immersed in hydrogen 
gas can exchange energy and particles with the environ- 
ment. In vacuum the exchange of energy is possible by 
contacting GO with a canonical thermostat that keeps 
temperature constant after reaching the thermodynamic 
equilibrium. In this case the exchange of particles is pos- 
sible either if the functional groups in GO are in a highly 
unstable configuration or if they undergo a chemical re- 
action and form chemical complexes. If the molecules 
formed by such complexes are stable they can be de- 
tached from their host and diffuse thermally away from 
GO. In an environment rich with hydrogen molecules we 
further expect that the reaction-diffusion processes are 
accounted for the particle exchange between GO and its 
environment. 

Our ab-initio/DFT calculation consists of Car- 
Parrinello molecular dynamics (CPMD) model^^, in 
which the potential energy of the system can be cal- 
culated on the fiy, as needed for the conformations of 
the dynamical trajectory to simulate the chemical reac- 
tion pathways. It is implemented in a plane- wave ba- 
sis within local spin density approximation (LSDA) with 
an energy cutoff of 70 Ry (rydberg), and with Becke^^ 
exchange and Lee- Yang-Parr (BLYP) gradient-corrected 
functional^^. Norm conserving ultrasoft Vanderbilt pseu- 
dopotentials were used for oxygen, hydrogen, and carbon. 
The CPMD microcanonical dynamics (constant energy 
ensemble) that allows achieving dynamical equilibration 
for T =300-1500 K were performed, after wave function 
optimization. Then a requenching of the wave function 
and canonical CPMD is performed. A supercell with pe- 
riodic boundary condition in the plane of graphene with 
Poisson solver of Martyna and Tuckerman was used^. 

Our approach for reactive MD consists of ReaxFI^ 
implemented in Large-scale Atomic/Molecular Massively 
Parallel Simulator (LAMMPS^P. The consistency and 
accuracy of ReaxFF with DFT calculation is checked sys- 
tematically by employing ab-initio CPMEB^. In ReaxFF 
the atomic interactions are described by the reactive force 
field potential^^. ReaxFF is a general bond-order de- 
pendent potential that provides accurate descriptions of 
bond breaking and bond formation. Recent simulations 
on a number of hydrocarbon-oxygen systems^^, GO^^ 
and organic molecules^^ showed that ReaxFF reliably 
provides energies, transition states, reaction pathways 
and reactivity trends in agreement with ab-initio cal- 
culations and experiments. The information on the 
type of chemical reactions and the time-evolution of the 




^0 100 200 300 400 500 600 
time (fs) 



FIG. 1: Epoxy and carbonyl pair formation of GO at 
T = 300K (a) and T = 1500K (b). (c) Kohn-Sham energy 
and C-C bond length (in carbonyl groups, shown by a dash- 
line in figure b) as a function of time and temperature. The 
transition state for C-C bond breaking corresponding to the 
energy barrier A ^ 3eV shown between the first two dash 
lines. A drop in Kohn-Sham energy is an indication of for- 
mation of carbonyl groups. The thermal fiuctuations in C-C 
bond length SA predicted by ab-initio CPMD shows the 
stability of carbonyl pair. The lowest energy difference be- 
tween the configurations in 1500K and 300K corresponding to 
pair of carbonyls and epoxys is calculated 0.09 eV within the 
present MD time. This small energy difference (almost one 
order of magnitude smaller than the one reported in Ref.!^ 
is due to relaxation of the structural strain by the surface 
diffusion of epoxys. Carbon and oxygen atoms are shown as 
green and red. 



damage is collected via running the MD up to 5 ps 
where the rearrangement of the atomic coordinates have 
been deduced from a dynamical trajectory calculated 
by ReaxFF. These simulations performed using periodic 
boundary conditions in a canonical NPT ensemble with 
a Nose-Hoover thermostat for temperature control and 
a time step of 0.25 fs. Finally a fitting to non-linear 
reaction-diffusion process is performed. 
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III. RESULTS AND DISCUSSIONS 



A. Carbonyl vs. epoxy formation, DFT results 



First we present series of ab-initio CPMD calculation 
to study the oxygen-graphene interaction as a function of 
temperature. At room temperature T = 300K, Fig. [TJa) 
shows formation of epoxy groups by two oxygen initially 
located in the opposite sides of the graphene sheet where 
carbon and oxygen atoms are shown with green and red 
balls. Increasing temperature up to T = 1500K, as shown 
in Fig. [TJb), change the pathway to carbonyl formation 
that accompany with formation of holes. Figs. [TJc) and 
(d) show the respective Kohn-Sham (KS) energies equiv- 
alent to potential energies in classical MD and carbon- 
carbon bond length (shown by arrow in Figs. [TJc) and 
(d)). The simulation is performed up to t = 600 fs where 
it shows stability of the system under study. Accordingly 
the transition for breaking C-C bond requires passing 
through an energy barrier A ^ 3eV. In high tempera- 
tures the shear force induced by anti-parallel and out-of- 
plane vibrational modes of oxygen-carbon bonds acting 
on carbon-carbon bonds provide enough kinetic energy 
that allows passing through the transition state. In this 
case the in-plane C-sp^ bond breaks and the formation of 
carbonyl pair lowers the local strain energy of graphene. 

In lower temperatures, however, the vibrational energy 
is not large enough to break the C-C bonds. Instead, sys- 
tem finds another pathway to a stable configuration by 
moving two oxygen one lattice constant away from each 
other such that it lowers the planar strain energy. This 
relaxation of local energy is assisted by the surface dif- 
fusion of the oxygen pairs. Because the chemical bond 
breaking is much faster process than mechanical expan- 
sion of the graphene lattice structure, the structural re- 
laxation in achieving the minimum energy configuration 
for a carbonyl pair occurs faster than epoxy pair. Within 
0.6 ps of CPMD time, we find the energy difference be- 
tween two lowest-energy configurations of carbonyl and 
epoxy pairs A^ks ^ —0.1 eV. However, a small and 
steady negative-slope seen in KS energy of epoxy pairs 
implies that the energy of epoxy pair asymptotically falls 
to a lower energy than the energy of carbonyl- pair. 

As the oxygenation of the graphene involves with 
passing through a large energy barrier in the reaction- 
diffusion processes, and because the ground state energy 
of both carbonyl and epoxy configurations are approxi- 
mately degenerate, the carbonyl formation is expected to 
be dominant at high temperature where the vibrational 
modes of the oxygen-carbon bonds fluctuate strongly and 
thus oxygen becomes more damaging whereas in lower 
temperature the dominant surface diffusion allows sta- 
bility of the epoxy groups without forming the structural 
holes. In the latter case the local strain energy can be 
relaxed through the thermal expansion of the lattice as- 
sisted by in-plane phonon modes. 
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FIG. 2: (a) Carbonyl formation at T = 1500K using 
ReaxFF-MD. (b) Similar simulation in H2 rich envoronment. 
The final products are water molecule and epoxy group, (c) 
C-C bond length (shown by dash line in figures a and b) as 
a function of time calculated at T = 1500K in vacuum and 
in H2 environment. Carbon, oxygen and hydrogen atoms are 
shown as green, red and white, respectively. 



B. Carbonyl formation, ReaxFF results 

The above analysis, consistent with the empirical 
data^^, indicates that formation of the carbonyl holes in 
high temperature annealing processes can not be avoided. 
A post-processing in a chemically active environment 
that allows removal of oxygen is an alternative possibil- 
ity that may lower the chance in carbonyl formation as 
well as reducing the oxygen concentration in GO. In this 
section, we examine such possibility. 

In Fig. [2|a) we describe a pathway in carbonyl pair 
formation in vacuum obtained after Ips MD simulation 
by ReaxFF. To check the consistency of the results with 
ab-initio CPMD, presented in previous section, we start 
with an initial configuration that consists of two oxygen 
located at the opposite side of graphene with separation 
of one lattice constant. Fig.[2|b) shows result of a similar 
MD simulation starting from the same initial configura- 
tion used for Fig. [2|a), except that hydrogen molecules 
are added in the computational box. Here carbon, oxy- 
gen and hydrogen atoms are shown as green, red and 
white balls respectively. 

We first energy minimize the system of hydrogen and 
leave a separation gap between GO and hydrogen gas to 
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FIG. 3: (a) An intermediate configuration corresponding to 
the onset of carbonyl formation found in the hydrogen-rich en- 
vironment. This configuration shows stability up to t ~ 500 
fs using ReaxFF-MD. Here two oxygen are initially located 
in the opposite sides of the graphene sheet and the temper- 
ature is T = 1500 K. (b) Final state of the same ReaxFF- 
MD that shows formation of hydroxyl and epoxy functional 
groups for t > 500 fs. This chemical pathway shows suppres- 
sion of carbonyl pair in favor of formation of hydroxyl and 
epoxy functional groups in the hydrogen-rich environment at 
t ~ 500 fs. Two dash-lines used to show the broken carbon- 
carbon bond length. (Bottom) C-C bond length as a function 
of time. The spike at the beginning is an indication of the on- 
set of carbonyl bond formation that is suppressed due to the 
presence of hydrogen molecules (seen as white spheres) and 
the transformation to hydroxyl and epoxy functional groups 
with dec ~ 1.6 A. Shown in inset is the same calculation in 
vacuum performed up to t = 5 ps. It shows a bi-st ability 
behavior in the equilibrium position of carbonyl-pair distance 
that jumps between 1.8 A and 2.5 A. 



avoid any bad contact and check the convergence of the 
results as a function of gap size to minimize the compu- 
tational artifact introduced by the separation gap. As 
it is shown, a water molecule forms after removal of one 
oxygen by hydrogen (for pictorial clarity, other hydrogen 
molecules are not shown). The second oxygen converts 
to an epoxy, hence the hole associated to carbonyl pair 
does not form. The bold-line in Fig. [2jc) shows the cor- 
responding carbon-carbon distance, dcc^ as a function of 
time in vacuum. At t = it starts from graphene lattice 
constant dec = 1.42 A and approximately at t = 50 fs it 
undergoes strong fluctuations and jumps to dec ~ 2.bA. 
In Fig. [2|c), the dash-line shows the time-evolution of 
C-C distance in the hydrogen-rich environment. We ob- 
serve the onset of carbonyl formation at the beginning in- 
dicated by a peak in the amplitude of C-C distance oscil- 
lations, however it is immediately suppressed by oxygen- 
hydrogen chemical interaction that leads to water forma- 
tion. 

Similar calculation based on ReaxFF-MD starting 
from slightly different initial condition is shown in Fig. [3] 
It illustrates an intermediate configuration and formation 
of carbonyl pair (Fig. [3|i) within t < 300 fs in which an 
absorption of hydrogen by oxygen occurs. This process 
converts one of the carbonyl pair oxygens to hydroxyl. 
Around t ~ 600 fs the second oxygen forms a bond with 
a neighboring carbon and an epoxy functional group is 
formed. This is shown in Fig. jsj). In Figjsjc), the time 
evolution of Cq-Cq bond length, corresponding to the 
dash line in the right side of Fig. [sja), shows gradual 
contraction, a manifestation on carbonyl pair suppres- 
sion. 

To check the establishment of carbonyl pair in vacuum 
for this particular simulation, we performed ReaxFF-MD 
up to 5 ps where we observed a bi-stability behavior in 
the equilibrium position of Cq-Cq distance that jumps 
between l.sA and 2.5A, seen in the inset of Fig [sj^c). 
Note that dec ^ 3A is predicted by ab-initio CPMD 
(see Fig. [ijc)). In vacuum, we observe that the onset 
of transition to Cq-Cq distance with 2.5A (carbonyl for- 
mation) occurs around t ~ 100 fs. This is in contrast 
with the same calculation in the hydrogen-rich environ- 
ment, in which Cq-Cq distance gradually decays to 1.6A 
(without forming carbonyl pair) until a hydroxyl func- 
tional group forms around t ~ 300fs. This comparison 
shows that within 100 < t < 300 fs, hydrogen-molecules 
in atmosphere hinder the carbonyl formation, by apply- 
ing a physical pressure to epoxy pair without forming a 
chemical reaction. 

We now turn to a system with large number of epoxy 
and hydroxyl functional groups intended to be a model 
for realistic structure of GO. Following Bagri et a/.^^^^, 
we start from an initial GO where the position of car- 
bons follow the long-range order structure of pristine 
graphene. Fig. [4] shows the initial and final structures 
in vacuum and in hydrogen-rich environment. The final 
structures are shown after Ips MD run. For clear visual- 
ization of the carbon structures, we omitted oxygens and 
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FIG. 4: Initial (left) and final states of GO after annealing 
at t = 2ps, and T = 1500K in vacuum (middle) and in H2- 
rich environment with the density of nn ~ 0.24 x 10^*^ cm"*^ 
(right). For clarity in visualization only carbons are shown 
in the final states. The holes in GO are due to formation 
carbonyl and other functional groups. 



hydrogens in the final configurations. We initially dis- 
tribute oxygens and hydrogens randomly as epoxy and 
hydroxyl groups on two sides of graphene sheet with ap- 
proximate area coverage of 12% epoxy and 7% hydroxyl 
using periodic boundary condition. Here the GO con- 
sists of Nc = 640, A^o = 122, and A^oh = 43, where 79 
oxygens are in form of epoxy and the rest are hydroxyl 
groups. 

The number of C-sp^ bonds per area is in particular the 
main factor that determines the quality of the graphene 
structure. These defects observed in GO, are absent in 
pristine graphene where such planar distortion does not 
occur and carbon coordination number does not change 
during MD by Hydrogen gas in environment. 

The number of carbons and oxygens bounded in 
graphene-sheet can be calculated as a function of time 



FIG. 5: Time evolution of the average coordinate number of 
carbon and oxygen calculated as average bond between two 
carbons (C-C) and carbon-oxygen (C-0). A cut-off distance 
= 2A is used to count the number of separated atoms from 
GO. 



by counting the number of C-C and C-0 bonds per car- 
bon and oxygen atom respectively. First we note that 
we found no chemically-unbounded carbon from the pla- 
nar GO either by thermal annealing or by hydrogen 
gas. Therefore we discard the possibility in the chem- 
ical removal of carbon from GO. The average number of 
carbon-bonds can be considered as representation of the 
carbon coordination number in the graphene-sheet. Dur- 
ing the annealing process, the epoxy and hydroxyl groups 
initially localized in the GO may interact and form free 
water molecules, OH free radicals, oxygen, ozone and hy- 
dronium molecules {H^O). To count the number of ab- 
straced oxygens, we calculate the distance between oxy- 
gens and carbons, doc- We define a length-scale as a 
cut-off such that if doc >d*=2A the chemical bonds is 
accounted for the broken-bonds. The result of this calcu- 
lation shown in Fig. |5] indicates that in vacuum the aver- 
age C-0 bonds reach steady state shortly after starting 
MD, however, in hydrogen-rich environment it decreases 
gradually as a function of time. We therefore conclude 
that hydrogen in environment do not suppress C-C bond 
breaking (e.g., etch hole formation), but that they fa- 
cilitate continuous removal of oxygen leading to reduced 
C-0 bonding in GO. This is consistent with the expec- 
tation on the role of hydrogen in GO reduction. 

Some interesting events responsible for the structural 
improvement of graphene sheet and observed in our sim- 
ulation include the collision between hydrogen molecules 
and oxygen functional groups that lead to capturing of 
hydrogen by oxygen and formation of a hydroxyl group 
on the surface of GO and finally formation of free wa- 
ter and hydronium molecules. In high temperatures 
where the average kinetic energy of hydrogen molecules 
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is larger, this happens more frequently and in a series 
of steps. For example, formation of transient hydrogen- 
oxygen complexes for a short period of time, accompanied 
by partial transfer of hydrogen kinetic energy to oxygen 
eventually lead to capturing of hydrogen by oxygen and 
forming OH on the surface of GO. 

Moreover, in regions of GO that the density of oxygen 
is high and they form cluster of oxygen, harmful chem- 
ical pathways by hydrogen can be dominant, in partic- 
ular in high temperatures. We can single out a process 
in which a hydrogen molecule dissociated by two neigh- 
boring oxygen atoms and form a pair of hydroxyl groups 
close to each other. Because this chemical reaction in- 
volves the energy transfer from hydrogen molecules to 
newly formed hydroxyl groups, the underlying C-C bonds 
undergo strong fluctuations that increases the chance of 
C-C bond breaking. With increasing the density of hy- 
drogen gas, oxygen localized on GO are nocked off with 
higher frequency allowing conversion of a pair of neigh- 
boring epoxys to carbonyl-holes and/or to hydroxyl- holes 
by collective C-C bond breaking. This implies an optimal 
limit on density of hydrogen gas to achieve the highest 
quality for graphene samples. 

We now turn to study the annealing process of GO 
that contains a large etch- hole. We further compare the 
annealing of GO with and without etch-hole. Figs. |6] 
shows the post-optimized initial structure of a single 
sheet of GO and the one with a large etch-hole. Un- 
like the previous calculation, the underlying structure 
of host-carbons is not a perfect periodic structure of 
graphene (because of structural optimization). For a 
comparison, we build a similar structure for GO (with- 
out etch-hole) and perform similar post-optimization be- 
fore starting thermal annealing. Here GO consists of 
Nc = 864, No = 255, and A/'oh = 144 corresponding 
to the number of carbon, oxygen (in form of mixture of 
carbonyl and epoxy groups) and hydroxyl-groups in GO, 
hence the oxygen percentage coverage is roughly 30%. 
For a GO with an etch- hole Nc = 836. To passivate 
the sp^ carbon-dangling bonds in the etch-hole, oxygen- 
atoms and hydroxyl groups are added hence the oxygen 
percentage coverage is slightly more than 30%. More- 
over, the average number of carbon-oxygen bonds per 
oxygen {Nq-o) ~ 1.6 and {Nq-o) ~ 1.5 for GO and 
GO with the etch-hole. This is because in the etch-hole 
the oxygen-carbon form double-bonds. The computa- 
tional box is filled up to 2000 hydrogen-atoms. Its planar 
and perpendicular dimensions are (4.3 x 5) x 2 nm^ with 
nn ^ 0.48 x lO^^cm"^ denoting the maximum density of 
hydrogen gas considered in our simulation. 

The time evolution of the oxygen concentration in GO 
as a function of temperature, measured by the average 
C-0 bonds, is shown in Fig. [7| As expected, the number 
of reduced oxygen increases with temperature. Compar- 
ing the time evolution of oxygen removed from GO and 
GO with a large etch-hole shown in Fig. [71(a) and (b), 
reveals that the initial oxygen abstraction close to large 
etch-holes take place with slower rate and without sig- 




FIG. 6: Initial structure of relaxed GO (top) and GO with 
a large etch-hole (bottom) used to study for the modeling of 
thermal annealing in the hydrogen-rich environment. These 
structures were first relaxed to their equilibrium structure 
then used for MD calculation. 
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FIG. 7: The time evolution of average carbon-oxygon bonds 
as a function of temperature and initial hydrogen density of 
riH ^ 0.24 X lO^^cm"^ for (a) GO and (b) GO with a single 
etch-hole. 

nificant dependence on temperature, most likely because 
oxygen near etch-hole are more strongly bounded to GO 
than oxygen in the bullP^. It is interesting to note that 
the dynamical behavior of the etch-hole interacting with 
hydrogen-gas is different in high temperatures, e.g., close 
to T = 1500 K the nucleation and growth of damaged 
C-C bonds in form of microscopic-cracks propagating to- 
ward the bulk is emerged. MD trajectories obtained by 
simulations of the etch-holes in vacuum and in the pres- 
ence of hydrogen gas at T = 1500 K are presented in 
on-line supplementary materials. To investigate the sta- 
bility of etch-hole in lower temperatures, a series of simu- 
lation in hydrogen gas at T = 1200 K with concentration 
up to flu ^ 0.8 X lO^^cm"^ (up to two times larger than 
the corresponding hydrogen gas concentration used at 
T = 1500 K) have been performed and no major damage 
in C-C bonds (as found at T = 1500 K) observed. 

The time-evolution of H-0 bonds is governed by non- 
linear reaction-diffusion process. In a chemical kinetics 
model, mass-action equations describe the interaction of 
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hydrogen molecules with GO to form complex, C*. This 
step is followed by a second-step of formation of a reac- 
tion product, P, and the release of new form of GO that 
contains less oxygen functional group 



[GO] + [H] ^ [C*] ^ [GO'] + [P], 



(1) 



with the reaction rates ki and /c2 associated with the first 
and second steps of chemical reactions. In this model we 
ignore the reverse reaction of C* GO + H2, assuming 
k-i « /c2. Also we assume that the kinetic model is 
reaction dominant such that the diffusion of hydrogen 
molecules occurs much faster than chemical reactions, 
hence the time-evolution of reaction-diffusion process is 
dominated by reaction rates. Denoting uqq and nco' 
as the initial and final concentration of oxygen in GO, 
and nc* and np as the concentration of all complexes 
and products, the mass-action equations corresponding 
to the reaction of equation [l] are given by 



dn 



GO 



dt 

dnc* 
dt 



dnco 
^ — ; — 
dt 



kinuTiGO - k2nc*, 



(2) 



subjected to the mass conservation nn -\-nc* = ^0 where 
no is the concentration of hydrogen molecules at t = 0. 
Here a is the oxygenation rate of graphene that lead for- 
mation of epoxy and hydroxyl functional groups on GO. 
To calculate solutions of Eqs. [l][2] numerically we intro- 
duce scaled variables g{t) = nQo{t) /uq^ c{t) = nc*{t)/no 
and ri = kiriQ. Using these new variables, Eqs. |2] trans- 
form to 



dc 



-c{tMt), 

c{t)]g{t) - k2c{t). 



(3) 



The numerical solutions of Eqs. [I]|3] provide a fitting 
model to our MD simulation shown in Figs. 7]|8 Fig. |9] 
shows the percentage oxygen remaining on GO and the 
fitting to the mass- action equations with ki = 0.15 and 
k2 = 0.1 ps-^ after annealing (a = 0) at T = 1200K for 
GO (with no hole). 



IV. CONCLUSION 

Ab-initio CPMD and ReaxFF molecular dynamic cal- 
culations performed in this work reveal that annealing 
of graphene-oxide in low-temperatures allows prevention 
of the structural damages due to relaxation of the lat- 
tice structure assisted by the oxygen surface diffusion. 
Our analysis indicates that the migration of epoxy func- 
tional groups provide a protection mechanism for C-sp^ 
bonds hence preserving the integrity of the graphene hon- 
eycomb structure. We further demonstrate that a chem- 
ically active environment, e.g. enriched with hydrogen 
gas, can assist removal of oxygen and blocks formation 
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FIG. 8: The annealing at T — 1200K as a function of density 
of H2 molecules shown in the figure in unit of lO^'^cm"^ for 
(a) GO and (b) GO with a single etch-hole. 
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FIG. 9: The percentage oxygen remaining on GO after an- 
nealing at T = 1200K as a function of density of H2 molecules 
shown in the figure in unit of lO^^cm"*^ for a GO with no 
etch-hole. The bold-lines are the result of fitting of MD (the 
dashed-lines with open circles and triangles) to the mass- 
action equations ([2|-([3| with a = 0, ki = 0.15 and fo 0.1 
ps"^ 



of carbonyl pairs, hence hindering the structural dam- 
age of GO either because of physical pressure exerted by 
the gas or imposed by chemical reactions. Our analysis 
shows that hydrogen in environment do not suppress G-G 
bond breaking, but that H facilitate continuous removal 
of oxygen leading to reduced G-0 bonding in GO, consis- 
tent with the expectation on the role of hydrogen in GO 
reduction. Finally, the results obtained in our studies 
can be used as a guideline for identifying the experimen- 
tal conditions that allows minimization of the mechanical 
and chemical damages to graphene and optimization of 
their performance. 
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